Exact quantum jump approach to open systems in Bosonic and spin baths 
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A general method is developed which enables the exact treatment of the non-Markovian quantum 
dynamics of open systems through a Monte Carlo simulation technique. The method is based on a 
stochastic formulation of the von Neumann equation of the composite system and employs a pair 
of product states following a Markovian random jump process. The performance of the method is 
illustrated by means of stochastic simulations of the dynamics of open systems interacting with a 
Bosonic reservoir at zero temperature and with a spin bath in the strong coupling regime. 
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The dynamics of open systems plays an important role 
in a wide variety of applications of quantum mechan- 
ics. The central goal of the physical theory |l| is to de- 
rive tractable equations of motion for the reduced density 
matrix ps{t) which is defined by the trace over the en- 
vironmental variables coupled to the open system. On 
the ground of the weak-coupling assumption the dynam- 
ics may be formulated in terms of a quantum dynamical 
semigroup which yields a Markovian master equation in 
Lindblad form However, strong couplings or inter- 
actions with low-temperature reservoirs give rise to large 
system-environment correlations which generally result 
in a failure of the Markov approximation: The open sys- 
tem dynamics develops long memory times and exhibits 
a pronounced non-Markovian behaviour. Important ex- 
amples of non-Markovian quantum phenomena are pro- 
vided by the atom laser [3| , the non-perturbative decay of 
atoms coupled to structured reservoirs environment- 
induced decoherence at low temperatures, and systems 
interacting with a spin bath 

The non-Markovian time-development of open quan- 
tum systems represents a challenge for the analytical 
and the numerical treatment, since it requires to cope 
with a complicated integro-differential equation 6], or 
with a highly non-local influence functional |2j. In the 
Markovian regime, Monte Carlo wave function tech- 
niques have been shown to provide efficient numerical 
tools 8J. In these techniques one propagates an en- 
semble of stochastic state vectors if){t) in the open sys- 
tem's state space such that the reduced density matrix is 
recovered through the ensemble average or expectation 
value ps(t) = E(\ip(t)){ip(t)\). In view of the numeri- 
cal efficiency of the method a generalization to the non- 
Markovian regime is highly desirable. One such general- 
ization || is based on a complicated stochastic integro- 
differential equation, involving a time-integration over 
the past history of the system dynamics. As shown in 
Ref. 0] one can av oid the solution of non-local equa- 
tions of motion by propagating a pair ipi{t), ip2(t) of 
stochastic state vectors and by representing the reduced 
system's density matrix through the expectation value 



ps(t) — E(\ijJi(t))(ip2(t)\). A similar idea has been used 
recently to design exact diffusion processes for systems 
of identical Bosons ^lj and Fermions [J^- However, 
this method demands that one first derives an appropri- 
ate time-local non-Markovian master equation for the re- 
duced system dynamics. Although this is indeed possible 
by use of the time-convolutionless (TCL) projection op- 
erator technique, the derivation becomes extremely com- 
plicated in the strong coupling regime which requires to 
take into account higher orders of the TCL expansion of 
the master equation. Another possibility [l3j is to exploit 
an explicit expression for the influence functional of the 
system. This method is, however, restricted to Gaussian 
reservoirs and linear dissipation. 

In this Letter an alternative approach to the open sys- 
tem dynamics is developed which is based on a stochastic 
formulation of the von Neumann equation for the density 
matrix p(t) of the total system. It is demonstrated that 
the motion of the state of any composite quantum sys- 
tem can be described through a Markovian random jump 
process which directly translates into a Monte Carlo sim- 
ulation technique for the full non-Markovian reduced sys- 
tem behaviour. By contrast to previous methods which 
mainly deal with Bosonic reservoirs and linear dissipa- 
tion, the present technique is generally applicable and 
does not require a specific form of the interaction be- 
tween system and environment, nor any assumption on 
the physical properties of the environment. 

In the interaction picture the Hamiltonian describing 
the system-environment coupling can be written as 



Hi{t) = ^A a (t) ® B a [t), 



(1) 



where A a (t) and B a {t) are system and environment op- 
erators, respectively. The dynamics of the total density 
matrix is then governed by the von Neumann equation, 

j t p(t) = -i[H I (t),p(t)]. (2) 
Our aim is to express p(t) through the mean value 

p(t) = E(|* 1 (i))($ a (i)|). (3) 
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|<E>i(i)) and |$2(f)) provide a pair of state vectors of the 
combined system which follow an appropriate stochastic 
time-evolution. These state vectors are taken to be direct 
products of system states ip v and environment states Xvi 



\® v {t))=il} v {t)®Xv{i), y=l,2. 



(4) 



Such a representation is indeed possible for any initial 
state yo(0) of the total system, regardless of the physical 
structure of open system and environment. In particular, 
p(0) needs not be an uncorrelated state. Equations © 
and Q amount to representing the density matrix of the 
open system through the expression 

Ps(t) = tr BP (t) = E (|Vi(*)><^(i)|<X2(t)|xi(*)» , (5) 

where tre denotes the trace over the variables of the envi- 
ronment. Contrary to the standard method, expression 
(J5J is an average over the product of \ipi)(4>2\ and the 
overlap (X2IX1) of the corresponding environment states. 

The task is to construct equations of motion for ij} v and 
Xv We suppose that the dynamics represents a piece- 
wise deterministic process (PDP). This is a Markovian 
stochastic process whose realizations consist of smooth 
deterministic pieces broken by instantaneous jumps oc- 
curring at random moments |l4j . A convenient way of 
formulating a PDP is to write stochastic differential equa- 
tions for the random variables. In our case an appropriate 
set of stochastic differential equations is given by 



dip v = ^2 



dxv = TyXvdt + 



A a - I j ij: v dN av , 

\xA\ 



\BaXv 

where / denotes the unit operator and 

\\A a ^ v \\ ■ \\B aXv 



(G) 



B a -I)xvdN av ,{7) 



I I'M I ■ WxA 



(8) 



Equations © and Q) relate the stochastic increments 
diji v = ip v {t + dt) - ij) v (t) and dxv = Xv{t + dt) - Xv(t) 
to a set of random integers dN av (t) which are known as 
Poisson increments and satisfy the relations 

dN av {t)dN^{t) = 5 a0 5 v ^dN av {t), (9) 
E(dN au {t)) = T av dt. (10) 

Equation © states that dN av (t) takes on the possible 
values or 1, while Eq. (fTT)f> tells us that dN av {t) = 1 
with probability T au dt. Moreover, under the condition 
that dN av (t) — 1 for a particular a and all other 
Poisson increments vanish. According to Eqs. 10 and (J7J) 
the state vectors then perform the instantaneous jumps 



"0. 



A a ipv, Xv 



\XA 



\B a Xv\ 



7B a Xv 



We observe that these jumps conserve the norm of the 
state vectors and occur at a rate T av defined in Eq. JSJ. 
Under the condition that all Poisson increments vanish, 
that is dN au (t) = for all a, v 1 we have dip v = and 
dxv = ^vXvdt. Consequently i/v remains unchanged dur- 
ing dt, while Xv follows a linear drift. Summarizing, ip^it) 
is a pure, norm-conserving jump process, whereas Xv(t) 
is a PDP with norm-conserving jumps and a linear drift. 

The stochastic differential equations © and 10 give 
rise to an equation of motion for the random matrix 
R(t) = |$i(i))($ 2 (i)l whose expectation value is re- 
quired to represent the total density matrix (see Eq. J3J)). 
To find this equation one applies the calculus of PDPs, 
which is analogous to the Ito calculus of the classical 
theory of Brownian motion. The calculations are partic- 
ularly simple in the present case since dN al/ dNp^ = 
for v 7^ /i (see Eq. (0). This implies that the in- 
crements and \d^2j are independent which yields 
dR = |rf$i)($ 2 | + \$i)(d§2\- Substituting the expres- 
sions (0J and using the stochastic differential equations 
© and 10 as well as the relations © and fTTHl one finds 

dR(t) = -i[H!(t),R(t)]dt + dS(t), (11) 
where dS = dT x R + RdT% is a stochastic increment and 
dT v = {? a vdt - dN av ) (I + iT-lA a B a ) . 

a 

Since the quantity T av dt — dN au is zero on average by 
virtue of Eq. l(TU)l. we conclude that E(dT u ) = and, 
hence, F,(dS) = 0. The role of the drift contributions 
T av dt is to compensate the mean contributions E(dN al/ ) 
of the jumps towards the stochastic increment dS. The 
average over the stochastic equation (jTT)l is therefore 
identical to the von Neumann equation (|2J. This shows 
that the PDP defined by Eqs. I© and indeed repro- 
duces the exact von Neumann dynamics of the combined 
system through the expectation value ©. 

The random process thus constructed leads to a Monte 
Carlo simulation technique of the open system dynamics. 
The numerical algorithms that can be used to generate 
a sample of realizations of the PDP are similar to the 
ones employed in the standard unraveling of quantum 
Markovian master equations. The decisive difference is, 
however, that the present method works with a pair of 
stochastic states and uses a representation in terms of 
states of the total system. Additionally, these features 
allow the determination of all kinds of multi-time quan- 
tum correlation functions directly from the Monte Carlo 
technique. 

An important limitation of the Monte Carlo method 
is set by the size a(t) of the statistical fluctuations. A 
detailed analysis of the stochastic differential equations 
reveals that a(t) is bounded for any finite time t, but 
may eventually grow exponentially at a rate of at most 
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2ro, where Tq is an upper bound for the transitions rates 
IV The stochastic simulation technique is thus useful for 
short and intermediate times, the relevant time scale be- 
ing given by l/2r . It is impossible, in general, to use 
the scheme for large times, trying to beat an exponential 
increase of the fluctuations by enlarging the number of 
realizations. However, the statistical errors may be con- 
siderably reduced by employing the fact that the |$„) 
evolve independently, or by using a class of stochastic 
states with a more complicated structure. 

To illustrate the method we first study the model of 
a two-state system with excited state |e), ground state 
\g), and transition frequency u>o- This system is coupled 
to a Bosonic reservoir consisting of a continuum of field 
modes k with creation and annihilation operators bt, 
and corresponding frequency ujk ■ The interaction Hamil- 
tonian is taken to be 

Hi(t) = <J+B(t) +a_5 t (i), (12) 

with the system operators a+ — \e)(g\ and er_ = \g)(e\ 
and the reservoir operator B(t) = Y^k 9kbke 1 ^ ^ u ' k ' t . 
The <7fe are mode-dependent coupling constants. We in- 
vestigate the dynamics evolving from the initial state 
1^(0)) = ip v (0) ® x*(0) = |e) <g> |0), where |0) denotes 
the vacuum state of the reservoir. The central quantity 
that determines the influence of the reservoir modes on 
the reduced system dynamics is provided by the correla- 
tion function (0\B(t')B^(t)\0) which may be expressed in 
terms of the reservoir spectral density J{oj). 



Q. 
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FIG. 1: The excited state probability p(t) of a two-state 
system coupled to a Bosonic reservoir: Analytical solution 
(continuous line) and Monte Carlo simulation of Eqs. @ , Q 
(diamonds) for 70/A = 5 with 2 ■ 10 7 realizations. The broken 
line shows the Born-Markov approximation. 

Figure n shows results of a Monte Carlo simulation of 
the stochastic differential equations © and Q for the 



damped Jaynes-Cummings model defined by the spec- 
tral density J(uS) = -^7 A 2 /[(w - ^) 2 + A 2 ]. The figure 
displays the population p(t) = (e\ps(t)\e) of the excited 
state, estimated with the help of expression (J5J. For the 
parameter values chosen the reservoir correlation time 
A -1 is five times larger than the reduced system's Marko- 
vian relaxation time 7^ . This shows that the case un- 
der study belongs to the strong coupling regime and that 
substantial deviations from the Born-Markov dynamics 
occur, as is indicated in Fig. ^ 

For small and intermediate couplings the reduced sys- 
tem dynamics is known to satisfy a non-Markovian mas- 
ter equation with a time-dependent generator, which can 
be obtained with the help of the TCL expansion in the 
coupling constant 70/A. However, the resulting pertur- 
bation expansion of the master equation diverges in the 
strong coupling regime 70/A > 0.5 for times which are 
larger than the first positive zero to of p(t). Beyond this 
singularity the generator of the master equation does not 
exist. The TCL master equation is therefore not capable 
of describing the revival of the excited state population 
for t > to. By contrast, the stochastic representation 
developed here reproduces the exact evolution of the re- 
duced system and, hence, correctly describes the dynam- 
ics even beyond the singularity. 

The stochastic simulation method can also be applied 
to the dynamics of open systems interacting with a spin 
bath, a quantum environment whose properties differ 
substantially from those of a Bosonic reservoir. We exam- 
ine a central spin model consisting of a spin with Pauli 
spin operator a which is coupled to a bath of N spins 
described by spin operators a&, j = 1,2,..., TV. The 
interaction Hamiltonian reads 

H x (t) =a 3 B 3 +a + B^(t) + a^B + {t), (13) 

where B±(t) = Y^j^^a ( £ e^ iuJot , B 3 = Ej^ 0)cr 3 J) - 
This model may be used to describe the interaction of 
a single electron spin, confined to a quantum dot, with 
an external magnetic field and a bath of nuclear spins 
through hyperfine interactions [l5j |. The transition fre- 
quency of the central spin is too and we set A® = A/VN, 
such that A provides the root mean square of the cou- 
plings of the central spin to the various bath spins. 

To give an example of a mixed initial state we choose 
p(Q) = \+)(—\®2~ N I, where |±) are eigenstates of the 3- 
component 03 of the central spin and / is the unit matrix 
in the state space of the bath. Initially, the bath is thus in 
a completely unpolarized state. This state can efficiently 
be realized through a mixture of basis states which are 
simultaneous eigenstates of the square of the total spin 
angular momentum J of the bath and of its 3-component 
J3, with an appropriate distribution of the correspond- 
ing quantum numbers j and m. Employing an argument 
given in 0] we conclude that the probability of finding 
the quantum numbers (j, m) in the initial mixture can be 
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written as P(j,m) = 2~ N ( ») - ( 



This 



N \ _ ( N 

KN/2+j] \N/2+j+lJ 

method of generating the initial bath state bears the ad- 
vantage that it enables one to employ the conservation 
of the 3-component ^(73 + J3 of the total spin angular 
momentum and to carry out a canonical transformation 
which removes the term 03-63 in the interaction Hamil- 
tonian (|13l) . This is an example of the implementation 
of known symmetries and conservation laws within the 
simulation algorithm. 
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FIG. 2: Real part 5Rp-i and imaginary part SJp^ of the 

coherence of a central spin coupled to a bath of 10 3 spins 
through the Hamiltonian (1131 with A/ujq = 0.5. Symbols: 
Simulation of Eqs. @, J7J with2-10 7 realizations. Continuous 
lines: Solution of the von Neumann equation JJ. Broken 
lines: TCL master equation in second order. 

Figure |21 shows the coherence p^ (t) = {+\ps(t)\— ) 

of the central spin obtained from a Monte Carlo sim- 
ulation employing again expression (0 for the reduced 
density matrix. To asses the performance of the stochas- 
tic method we compare the simulation results with the 
solution of the von Neumann equation for the total sys- 
tem. The simulation of the PDP obviously reproduces 
the von Neumann dynamics with high accuracy, the sta- 
tistical errors of the simulation being smaller than the 
size of the symbols. For the parameter values chosen, 
the dynamics significantly deviates from the predictions 
of the second order TCL master equation. 

The stochastic formulation naturally lends itself to sev- 
eral interesting generalizations. One possibility is to for- 
mulate the dynamics in imaginary time by means of the 
substitution t — > — i(3, to examine the equilibrium prop- 
erties of the system (/3 denotes the inverse temperature) . 
With slight modifications of the differential equations © 
and J7J) one finds a stochastic dynamics in imaginary time 
which unravels the canonical density matrix of the sys- 
tem via the expectation value p(/3) = E(|$i( / 9))($ 2 ( / 9)|). 

It was assumed in Eq. Q that the 1$^) are tensor prod- 
ucts of state vectors of system and environment. Since 



the interaction creates system-environment correlations 
it may be advantageous to employ a certain class of en- 
tangled stochastic states with the aim of a more efficient 
representation of p(t) as mean value over the underly- 
ing random process. A further possibility is to use a 
stochastic propagation of mixed states. This is achieved, 
for example, by the introduction of a stochastic matrix 
R(t) — |V'i(t))('02(t)| <8> -Rb(*)j where the ip v are states 
of the open system and i?B is a random operator in the 
state space of the environment. It is then again possible 
to construct stochastic differential equations leading to 
the exact von Neumann dynamics with the help of the 
mean value pit) = E(i?(i)). A systematic exploration 
of these ideas could be of great relevance for the devel- 
opment of efficient numerical algorithm of the quantum 
dynamics of open systems. 
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